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Abstract. We consider the problem of state and parameter estimation for a class of nonlinear 
oscillators defined as a system of coupled nonlinear ordinary differential equations. Observable 
variables are limited to a few components of state vector and an input signal. This class of systems 
describes a set of canonic models governing the dynamics of evoked potential in neural mem- 
branes, including Hodgkin-Huxley, Hindmarsh-Rose, FitzHugh-Nagumo, and Morris-Lecar mod- 
els. We consider the problem of state and parameter reconstruction for these models within the 
classical framework of observer design. This framework offers computationally-efficient solutions 
to the problem of state and parameter reconstruction of a system of nonlinear differential equa- 
tions, provided that these equations are in the so-called adaptive observer canonic form. We show 
that despite typical neural oscillators being locally observable they are not in the adaptive canonic 
observer form. Furthermore, we show that no parameter-independent diffeomorphism exists such 
that the original equations of these models can be transformed into the adaptive canonic observer 
form. We demonstrate, however, that for the class of Hindmarsh-Rose and FitzHugh-Nagumo 
models, parameter-dependent coordinate transformations can be used to render these systems into 
the adaptive observer canonical form. This allows reconstruction, at least partially and up to a 
(bi)linear transformation, of unknown state and parameter values with exponential rate of conver- 
gence. In order to avoid the problem of only partial reconstruction and at the same time to be able 
to deal with more general nonlinear models in which the unknown parameters enter the system 
nonlinearly, we present a new method for state and parameter reconstruction for these systems. 
The method combines advantages of standard Lyapunov-based design with more flexible design 
and analysis techniques based on the notions of positive invariance and small-gain theorems. We 



'Corresponding author. E-mail; I.Tyukin@le.ac.uk 



1 



D. Fairhurst et al. 



Observers for Canonic Models 



show that this flexibility allows to overcome ill-conditioning and non-uniqueness issues arising in 
this problem. Effectiveness of our method is illustrated with simple numerical examples. 

Key words: Parameter estimation, adaptive observers, nonlinear parametrization, convergence, 
nonlinear systems, neural oscillators 

AMS subject classification: 93B30, 93B10, 93B07, 93A30, 92B05 



Notations and Nomenclature 

The following notational conventions are used throughout the paper: 

• M is the field of real numbers. 

• M>o = {x e M I X > 0}. 

• Z denotes the set of integers, and N stands for the set of positive integers. 

• The Euclidian norm of x G M" is denoted by . 

• C denotes the space of continuous functions that are at least r times differentiable. 

• Let /i : — M be a differentiable function and / : M" — )• M". Then Lfh{x), or simply 
Lfh, is the Lie derivative of h with respect to /: 

• Let /, (jf : M" — )■ M" be differentiable vector- fields. Then the symbol [/, g] stands for the Lie 
bracket: 

ox ox 

The adjoint representation of the Lie bracket is defined as 

• Let ^ be a subset of M", then for all x E M", we define dist(^, x) = inf^g^ II — q\\ ■ 

• 0{-) denotes a function such that lim^^o / s = d, d E M, d ^ 0. 

• Finally, let e E ffi>o, then ||x||^ stands for the following: 



\x\ 



0, ||x|| < e. 
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1. Introduction 

Mathematical modelling of brain processes and function is recognized as an important tool of mod- 
em neuroscience [fT4|. It allows us to predict, analyze and understand intricate processes of neural 
computations, without invoking technically involving and costly experiments. Successful exam- 
ples include but are not limited to modelling the memory function [fTSl and the mechanisms of 
phase-resetting in olivo-cerebellar networks |[38l . Availability of quantitatively accurate models of 
individual neural cells is an important prerequisite of such studies. 

The majority of available models of individual biological neurons are the systems of ordinary 
differential equations describing the cell's response to stimulation; their parameters characterize 
variables such as time constants, conductances, and response thresholds, important for relating 
the model responses to behavior of biological cells. Typically two general classes of models 
co-exist: phenomenological and mathematical ones. Models of the first class, such as e.g. the 
Hodgkin-Huxley equations, claim biological plausibility, whereas models of the second class are 
more abstract mathematical reductions without explicit relation of all of their variables to physical 
quantities such as conductances and ionic currents (see Table [U). 

Despite these differences, these models admit a common general description which will be 
referred to as canonic. In particular, the dynamics of a typical neuron are governed by the following 



Table 1 : Examples of typical mathematical models of spiking single neurons. Biologically plau- 
sible equations of membrane potential generation in a giant axon of a squid (left panel), and a 
reduction of these equations to an oscillator with polynomial right-hand side (right panel) 



Hodgkin-Huxley Model [fTOl] 


Hindmarsh-Rose Model S 


V = I{t) - [d^m^h{v + 02) + e:in\v + 6^) + O^v + 6^] 
m = (1 — mjv? 1 — - — 1 — m 6*9 exp — 

V ^8 / ^10 

n = (1 n)Qx\^\ . ) n 6114 exp (1.1) 

h = (1 - h)eie exp-^-h/ ( 1 + exp ^ "j" ) 

ip{x) = x/{expx — 1), 9i, ...,6ig - parameters, J : M — )■ M - 
input current 


V = ei{e2r-f{v)) + i{t) 

r = e:,{g{v) - e^r) (1.2) 

f{v) and g{v) are polynomi- 
als: 

f{v) = e, + e,v + 

g{v) = eg + e^ov + OnV^ 

/ : M — !• R is the external in- 
put, 9i, ...,9u- parameters 
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set of equations 



V 



-ai{v,9,t)ri + bi{v,9,t); 
{01,02,...). 



(1.3) 







in which the variable v is the membrane potential, and are the gating variables of which the 
values are not available for direct observation. Functions (pj{., ■), pj{.) E are assumed to 
be known; they model components of specific ionic conductances. Functions aj(-), G 
are also known, yet they depend on the unknown parameter vector 9. System (11.31) is a typical 
conductance-based description of the evoked potential generation in neural membranes |fT3l . It is 
also an obvious generalization of many purely mathematical models of spike generation such as the 
FitzHugh-Nagumo |I3 or the Hindmarsh-Rose equations |i9j|. In this sense systems (ll.3|) represent 
typical building blocks in the modelling literature. 

In order to be able to model the behavior of large numbers of individual cells of which the 
input-output responses are described by (11.31 ). computational tools for automated fitting of models 
of neurons to data are needed. These tools are the algorithms for state and parameter reconstruction 
of ( 11.31) from the available measurements of v{t) and J(t) over time. 

Fitting parameters of nonlinear ordinary differential equations to data is recognized as a hard 
computational problem |I51 that "has not been yet treated in full generality" EOl . Within the field of 
neuroscience, conventional methods for fitting parameters of model neurons to measured data are 
often restricted to hand-tuning or exhaustive trial-and-error search in the space of model parameters 
||3T| . Even though these strategies allow careful and detailed exploration in the space of parameters 
they suffer from the same problem - the curse of dimensionality. 

Available alternatives, recognizing obvious nonlinearity of the original problem, propose to 
reformulate the original estimation problem as that of searching for the parameters of a system 
of difference equations approximating solutions of ( 11.31) [IJ; or predominantly offer search-based 
optimization heuristics (see [37] for a detailed review) as the main tool for automated fitting of 
neural models. Straightforward exhaustive-search approaches however are limited to varying only 
few model parameters over sparse grids, e.g. as in OTIl where 8 parameters were split into 6 
bands. Coarseness of this parametrization leads to non-uniqueness of signal representation, leaving 
room for uncertainty and inability to distinguish between subtle changes in the cell. More fine- 
grained search algorithms are currently infeasible, technically speaking. Other heuristics, such 
as evolutionary algorithms, are examined in According to O, replacing exhaustive search 
with evolutionary algorithms allows to increase the number of varying parameters to 24. Yet, 
computational complexity of the problem still delimits the search to sparse grids (6 bands per single 
parameter) and requires days of simulation by a cluster of 10 Apple 2.3 GHz nodes. Furthermore, 
because all these strategies are heuristic, accuracy of final results is not guaranteed. 

The main aim of this article is to present a feasible substitute to these heuristic strategies for 
automatic reconstruction of state and parameters of canonic neural models (11.31) . To develop com- 
putationally efficient procedures for state and parameter reconstruction of (11.31) we propose to 
exploit the wealth of system-identification and estimation approaches from the domain of control 
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theory. These approaches are based on the system-theoretic concepts of observability and identifi- 
ability [[301 . [[T2l . [[T9l from control theory, and the notions of Lyapunov stability [[221 and weakly 
attracting sets [[261 . The advantage of using these approaches is that there is an abundance of al- 
gorithms (observers) already developed within the domain of control. These algorithms guarantee 
asymptotic and stable reconstruction of unmeasured quantities from the available observations, 
provided that the system equations are in an adaptive observer canonical form. Moreover, this 
reconstruction can be made exponentially fast without the need of substantial computational re- 
courses. We study if system ([1.3[) is at all observable with respect to the output v, that is if its 
state and parameters can be reconstructed from observations of v. We present and analyze typi- 
cal algorithms (adaptive observers) that are available in the literature. We show that for a large 
class of mathematical models of neural oscillators at least a part of the model parameters can be 
reconstructed exponentially fast. 

In order to deal with more general classes of models and also to recover the rest of the model 
parameters we introduce a novel observer scheme. This scheme benefits from 1) the efficiency of 
uniformly converging estimation procedures (stable observers), 2) success of explorative search 
strategies in global optimization by allowing unstable convergence along dense trajectories, and 3) 
the power of qualitative analysis of dynamical systems. We present a general description of this 
observer and list its asymptotic properties. The theory of this new class of algorithms is based 
on the results of our previous studies in the domain of unstable convergence [[351 , [[361 . We will 
present examples to demonstrate the performance of these algorithms. 

The paper is organized as follows. In Section 2 we provide the basic notions of observability 
from the domain of mathematical control, test if typical canonical neural oscillators are observable, 
and present two major classes of systems (canonical forms) for which computationally efficient 
reconstruction procedures are available. In Section 3 we analyze the applicability of standard 
observers to the problem of reconstructing all unmeasured variables and parameters of typical 
models of neurons. We present two special cases in which such reconstruction is possible. In 
Section 4 we provide a description and asymptotic properties of our algorithm that applies to the 
most general subset of models ([1.3[ ). Section [5] contains examples of application of the considered 
observers, and Section 5 concludes the paper. Proofs of the main technical statements are presented 
in the Appendix. 

2. Observer-based approaches to the problem of state and pa- 
rameter estimation 

Let us consider the following class of dynamical systems 



X = f{x,e) + g{x,e)u{t), xito) eVL^dW 
y = h{x), X eW, e eR'^, y eR 



n 



(2.1) 
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where /, 51 : x R'" -> M", /i : M" ^ M are smooth function^, and m : R ^ M. Variable x 
stands for the state vector, m G W C [to, 00) is the known input, E is the vector of unknown 
parameters, and y is the output of (12.11) . System (12.11) includes equations (11.31 ) as a subclass and 
in this respect can be considered as plausible generalizations. Obviously, conclusions about (12.11) 
should be valid for systems (11.31) as well. 

Given that the right-hand side of (12.11) is differentiable, for any x' G ^Ix, u G C^[to, 00) there 
exists a time interval T = [to, ti], ti > to such that a solution x{t, x') of (12.11 ) passing through x' at 
to exists for all t E T. Hence y{t) = h{x{t)) is defined for all t E T. For the sake of convenience 
we will assume that the interval T of the solutions is large enough or even coincides with [to, 00) 
when necessary. 

We are interested in finding an answer to the following question: suppose that we are able 
to measure the values of y{t) and u{t) precisely; wether and how the values of x' and parameter 
vector 9 can be recovered from the observations of y{t) and u(t) over a finite subinterval of T? A 
natural framework to answer to these questions is offered by the concept of observability ll30l . 



Definition 1 (Observability). Two states xi,X2 E R" are said to be indistinguishable (denoted 
by X1IX2) for ( I2.il) if for every admissible input function u the output function t — )■ y{t, 0, xi, u), 
t > Q of the system for initial state x(0) = Xi, and the output function t — y(t, 0, X2, u), t > 
of the system for initial state x{0) = X2, are identical on their common domain of definition. The 
system is called observable if XiXx2 implies Xi = X2. 

According to Definition [H observability of a dynamical system implies that the values of its 
state, x(t), t G [ti,t2] are completely determined by inputs and outputs u{t), y(t) over [ti,t2]. 
Although this definition does not account for any unknown parameter vectors, one can easily see 
that the very same definition can be used for parameterized systems as well. Indeed, extending 
original equations (12.11) by including parameter vector 6* as a component of the extended state 
vector X = (x, 9)-^ results in 

X = f{x,9)+g{x,9)u{t), 

9 = (2.2) 



y = h{x), x(to) G fi^ C 



or, similarly, in 



X = fix) +g{x)u{t) 

y = h{x), x(to) = (a;(to), 9f E C R'^+'^ 

where f{x) = {f{x, 9), 0)^, g{x) = {g{x, 9), 0)^, and h{x) = {h{x), 0). All uncertainties in (12.11) . 
(12.21) . including the parameter vector 9, are now combined into the state vector of (12.31) . Hence 
the problem of state and parameter reconstruction of (12.11) can be viewed as that of recovering the 
values of state for (12.31 ). 



^Let us recall that a function is smooth in G if for every x e G and n e iV the function -^f{x) is always defined. 
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Definition \T\ characterizes observability as a global property of a dynamical system. Some- 
times, however, global observability of a system in is not necessarily needed. Instead of asking 
if every point in the system's state space is distinguishable from any other point it may be suf- 
ficient to know if the system's states are distinguishable in some neighborhood of a given point. 
This necessitates the notion of local observability [[30l . 

Let V be an open subset of M". Two states xi,X2 G V are said to be indistinguishable (denoted 
by XiZ^X2) on V for (12.11 ) if for every admissible input function n : [0, T] — M with the property 
that the solutions x{t, 0, Xi, u), and x{t^ 0, X2, u) both remain in V for t < T the output function 
t — )■ 0, xi, m), t > of the system for initial state x(0) = xi, and the output function t — )• 
y{t, 0, X2, m), t > of the system for initial state a;(0) = X2, are identical for < t < T on their 
common domain of definition. 

Definition 2 (Local observability [{30l). The system is called locally observable at xq if there exists 
a neighborhood W dW^ ofxo such that for every neighborhood V <ZW ofxo the relation xqI^xi 
implies Xq = Xi. The system is locally observable if it is observable at each Xq. 

A number of observability tests are available that, given the functions h, f in the right-hand side 
of (12.31) . indicate if a given system is observable. Particular formulations of these tests may vary 
depending on whether e.g. the functions /, g are analytic or time-invariant (inputs are constants). 

In this article we will restrict our attention to those systems (12.21 ) in which the inputs u(t) are 
constants. In this case we can replace the function u{t) with an unknown parameter, and system 
(12.21) can be viewed as a system (12.31) yet without inputs. One of the most common observability 
tests for this class of autonomous systems is given below (see also |[30ll . Theorem 3.32): 

Proposition 3 (Observability test (Corollary 3.33, Il30l )). System {\2.3\l is locally observable at a 
point x° eU C M"+'^ if 

rank^(h{S;) L^-h{x) Lyh{x) ... L^+'^-^/i(5))^ = n + VxGf/ (2.4) 

In what follows we shall use the test above to determine if the models of neural dynamics are 
at all observable. 

2.1. Local observability of neural oscillators 

We start our observability analysis by applying the local observability test (12.41) to the Hindmarsh- 
Rose model (11.21) . In order to do so we shall extend the system state space so that unknown 
parameters are the components of the extended state vector. In the case of the Hindmarsh-Rose 
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(a) 6 = 10^0 (b) 6 = IQis (c) S = 10^ 

Figure 1 : Observability tests for Hindmrsh-Rose model neuron (11.21) 



model this procedure leads to the following extended system of equations: 



6'l3 
6*12 

^22 
6*21 

V ^ / 



/ 9i3xl + Ouxl + 9nXi + 6^10 + X2 \ 
-Xx2 + O22XI + 621X1 










(2.5) 



To test if there are points of local observability of system (12.51) it is sufficient to find a point in 
the state space of (12.51) at which the rank condition (12.41) holds. Here we computed the determinant: 



d 

D{x) = — Lfh{S;) L)h{i) ... 

X = {xi, X2, 6^13, 612, Oil, ^10, ^22, ^21, 



T n—l 



h{i)y 



on a sparse grid (of 101 x 101 pixels) and plotted those regions for which the determinant is less 
than a certain value, 5. The neuron parameters were set to L = —2, Oi^ = —10, 6*12 = —4, On = 
6, 6'io = 1,6^22 = —32, 6^21 = —32. Figure [T] shows results (obtained using Maple) for various 
values of 5. The shaded regions correspond to the domains where D{x) < 5. According to these 
results, when the value of delta is made sufficiently small, condition D(x) > 5 holds for almost 
all points in the grid. This suggests that there are domains in which model (11.21) is indeed at least 
locally observable. 

Let us now consider a more realistic, with respect to biological plausibility, set of equations. 
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(a) S = 10000 (b) S = 1000 (c) 6 = 100 

Figure 2: Observability tests for Morris-Lecar model neuron 



One of the simplest models of this type is the Morris-Lecar system ||28l : 



v{t) = -9ca (l + itanh (^) - E 



El = 100, E2 



wit) 
-70,Es 



- 9Kw{t){v{t) 
itanh (^^'^-wit) 

-50, = 15,^5 = 30,^6 = 60, (?ca 



E2) - gmi<t) - Es) 



cosh ( ^ 



(2.6) 



1.1, gx = 2.0, gr, 



0.5 



As in the previous example we extend the system state space by considering unknown parameter 
as components of the extended state vector. This extension procedure results in the following set 
of equations: 



v 



v{t) 
wit) 

9Ca{t) 
QKit) 

gmit) 

m 



\ 



( -gca (I + I tanh (^) [v - 100)) - gkw{v + 70) - (7, 
|(i + itanh(^)-^) cosh(|j) 




' + 50) \ 



V 










For this extended set of equations we estimated the regions where value of D{x) exceeds some 
given 5 > 0. These regions for different values of 5 are presented in figure [2] These results demon- 
strate that the Morris-Lecar system (12.61) is also locally observable. 

As we have seen above, a fairly wide class of canonical mathematical and conductance-based 
models of evoked responses in neural membranes satisfy local observability conditions. We may 
thus expect to be able to solve the reconstruction problem for these models. In fact, as we show 
below in Sections 3, 4, the reconstruction problem can indeed be resolved efficiently at least for a 
part of unmeasured variables of the system. However, before we proceed with detailed description 
of these reconstruction algorithms, let us first review classes of systems for which solutions to the 
problem of exponentially fast reconstruction of all components of state and parameter vectors are 
already available in the literature. 
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nxp 



2.2. Bastin-Gevers canonical form 

We start with a class of systems comprising of a linear time-invariant part of which the equations 
are known and an additive time-varying component with linear parametrization. Parameters of this 
time-varying component are assumed to be uncertain. This class of systems was presented by G. 
Bastin and M. Gevers in 1989, [|3||, and its general form is as follows: 

x = Rx + n{t)e + git) 
y{t) = xi{t) 

In (12.71) . a; G M" is the state vector with y = xi assigned to be the output. 9 eW = {9i,- ■ 9pf 
is the vector of unknown parameters, i? is a known matrix of constants where k'^ = (/c2, ■ ■ ■, ^; 
and F has dimension (n — 1) x (n — 1) with eigenvalues in the open left half plane. G 
is an n X p matrix of known functions of t; the first row is designated fii and the remaining n — 1 
rows f2. The vector function g{t) : M — )■ M" is known. 

Equations (12.71) are often referred to as an adaptive observer canonical form. This is because, 
subject to some mild non-degeneracy conditions, it is always possible to reconstruct the vector of 
unknown parameters 9 and state x from observations of y over time. Moreover, the reconstruction 
can be made exponentially fast. Shown below is the adaptive observer presented in [[3l. The 
system to be observed, state estimator, parameter adaption, auxiliary filter and regressor are given 
in equations (Q, ^M, ^M, (l2?T0l) . (ETTl) respectively 

£{t) = Rx{t) + Q{t)9 + g{t)+(^^'^^^.^ (2.8) 

§{t) = r^{mt)_ (2.9) 

V(t) = FV{t) + n{t), 1/(0) =0 (2.10) 

cp{t) = v^{t)k + nj{t) (2.11) 

The output isy = xi, its estimate isy = xi and its error isy = y — y. This observer contains some 
parameters of its own which are at the design's disposal. F = is an arbitrary positive definite 
matrix, normally chosen as F = diag(7i, 72, 7p), 7^ > 0. Ci > 0. The auxiliary filter V{t) is an 
(n — 1) X p matrix and (^(t) is ap vector. 

Using the transformation (12.121) . the error system (12.131) is obtained. 

X* = X — ( y~ j , X = X — X (2.12) 



-ci F 
F 

9 = - 



X 



ip^9 




(2.13) 
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It is shown in [3], for constant unknown parameters, that the solution x{t) = x{t), 9 = 6(t) 
of the extended system ( 12.81 ). (12.91 ). (12.101) . ( 12.111) is globally exponentially stable provided certain 
conditions on the regressor vector, tp{t), are met. These conditions are: 

• the regressor vector ipit) is bounded for all t > 

• is bounded for alH > except possibly at a countable number of points {tj} such that 
min|tj — tj l > A > for some arbitrary fixed A. 

• Lp(t) is persistently exciting: that is, there exists positive constant a, T such that for all to > 

to+T 

ip{t)ip'^{t)dt> al > (2.14) 
Formally, asymptotic properties of observer (12.81) . (12.91 ) are specified in the theorem below ll3lj^ 



t, 



Theorem 4. Suppose that 

1) ci > and F is a Hurwitz matrix, that is its eigenvalues belong to the left half of the complex 
plane; 

2) the function ipit) is globally bounded in t, and its time derivative exists and is globally 
bounded for all t > 0; 

3) the function Lp{t) is persistently exciting. 

Then the origin of ( 12.731) is globally exponentially asymptotically stable. 

Adaptive observer canonical form (12.71) applies to systems in which the regressor VL{t)6 does 
not depend explicitly on the unmeasured components of the state vector. The question, however, 
is when a rather general nonlinear system can be transformed into the proposed canonical form. 
This question was addressed in |[24l in which a modified adaptive observer canonical form was 
proposed together with necessary and sufficient conditions describing when a given system can be 
transformed into such form via a diffeomorphic coordinate transformation. This canonical form is 
described in the next subsection. 

^Here we provide a slightly reduced formulation of the main statement of ||3] corresponding to the case in which 
the values of 6 do not change over time. 



11 



D. Fairhurst et al. 



Observers for Canonic Models 



2.3. Marino-Tomei canonical form 



The canonical form presented in [24] is now shown here. The system to be observed (12.151) . state 
estimator (12.161) and parameter adaption (12.171 ) are given below 

x{t) = A^x{t) + Myit),u{t)) + bJ2^^^,/3i{y{t),u{t))ei 
Cix{t) 



y{t) 



1 ° 


1 


. 


.. 0\ 








1 . 


.. 








.. 








. 


.. 1 







. 


.. y 







.. 


0) 



(2.15) 



X 



GM",2/(t) GM,A(-,-) 

In (|2.15l) x{t) G is the state vector with xi{t) assigned to be the output y. Matrices Ai, Ci are 
in canonical observer form. 9 G = (6'i, ■ ■ 9p)^ is the vector of unknown parameters. The 
functions /3j are known, bounded and piecewise continuous functions of y{t), u{t). The column 
vector 6 G is assumed to be HurwitzQ with 6i 7^ 0. 
Shown below is the adaptive observer presented in Il24ll 



e 

K 



(Al - KCi)x{t) + 0o(l/(t), M(t)) + h J2 Hy{t),u{t))e, + Kyit) (2.16) 



i=l 



r/3(t)(y-Cix)sign(6i) 
^{Aih + \h) = {h,---,knf 



(2.17) 
(2.18) 



with r an arbitrary symmetric positive definite matrix and A an arbitrary positive real. The n x 1 
vector, h, is Hurwitz with 6„ 7^ 0. 

For the more general case where the vector h is an arbitrary vector, an observer is presented in 

m. 

Theorem 5. / |2?1/ There exists a local change of coordinates, z = $(x), transforming 



X 



x) 



y = xi 



4 = 1 



X eW\y eR,6i e M, qi : M" 
with h{x°) = and (/, h) an observable pair, into the system 



,n>2 



(2.19) 



i=l 



z G 



(2.20) 



"^We say that a vector b = (61, . . . , 6„)-'^ S M" is Hurwitz if all roots of the corresponding polynomial bip 
+ bn^ip + bn have negative real part. 



,n-l 
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with {Al, Ci) in canonical observer form d2.i5P . if and only if 

(i) [ad'jg, acf^g] = 0, < i, j < n — 1 

(ii) [qi, ad}g] = 0, < j < n - 2, l<i<p 
where the vector field, g{x), is uniquely defined by 



d_ 

dx 





h{x) 






1 \ 




Lfh{x) 

















\ 


LY^h{x) 


I 




1 1 / 



(2.21) 



The proof of this result is made along the following lines. Suppose we use the change of 
coordinates: z = then we have 



dx 



X 



9$ 

—f{x, 

9x dx 



(2.22) 
(2.23) 



It is shown in ||23l that providing we meet the constraint: 

[ad'jrg, ad^j^g] = 0, Wx e U,0 < i, j < n - 1 
then we can cast the system into the adaptive observer canonical form 

z = Aiz + 2p{y) + 2_^—9iqi{x), y = Ciz 



(2.24) 



i=l 



(2.25) 



with Al, Ci in canonical observer form (12.151) . Furthermore, it is shown in ll23l that providing we 
meet the constraint 



[qi, ad^^g] = 0, \/x e U,0 < j < n - 2 
then we can put the system into 



(2.26) 



Aiz + ^{y) + J2^My), y = Ciz 



(2.27) 



i=l 



This representation is linear in the unknown variables, Zi{t), Z2{t), Zn{t), and 9i{t),92{t), ...,9p{t), 
while it is nonlinear only in the output, y{t), which is available for measurement. 
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3. Feasibility of conventional adaptive observer canonical forms 

In this section we consider technical difficulties preventing straightforward application of conven- 
tional adaptive observers for solving the state and parameter reconstruction problems for typical 
neural oscillators. We start with the most simple polynomial systems such as the Hindmarsh-Rose 
equations. We show that even for this relatively simple class of linearly parameterized models the 
problem of reconstructing all parameters of the system is a difficult theoretical challenge. Whether 
complete reconstruction is possible depends substantially on what part of the system's right-hand 
side is corrupted with uncertainties. Despite in the most general case reconstruction of all compo- 
nents of the parameter vector by using standard techniques may not be possible, in some special 
yet relevant cases estimation of a part of the model parameters is still achievable in principle. 

Let us consider, for example, the problem of fitting parameters of the conventional Hindmarsh- 
Rose oscillator to measured data. In particular we wish to be able to model a single spike from 
the measured train of spikes evoked by a constant current injection. Classical two-dimensional 
Hindmarsh-Rose model is defined by the following system: 



in which / G M stands for the stimulation current. Trajectories x{t) of this model are known to 
be able to reproduce a wide range of typical responses of actual neurons qualitatively. Quantita- 
tive modelling, however, requires the availability of a linear transformation of (x{t),y(t)) so the 
amplitude and the frequency of oscillations x{t) can be made consistent with data. 

In what follows we will consider (13.11 ) subject to the following class of transformations: 



where ki > and Px, Py E M are unknown. Transformations (13.21) include stretching and trans- 
lations as a special case. In addition to (13.21) we will also allow that the time constants in the 
right-hand side of (|3.1I) be slowly time-varying. This will allow us to adjust scaling of the system 
trajectories with respect to time. 

Taking these considerations into account we obtain the following re-parameterized description 



X = —ax^ + hx^ + y + I 

y = c~ dx^ — y, a = 1, b = 3, c = 1, d = 5, 



(3.1) 




(3.2) 



of model (1311): 



3 




i=0 



(3.3) 



3 




1=1 



Alternatively, in vector-matrix notation we obtain: 




(3.4) 



14 



D. Fairhurst et al. 



Observers for Canonic Models 



where 

, , , , / 1 \ , , / xi 1 \ 

= ^ -A J ' = ^ xl xl xj (3.5) 

^ = (^1,3; ^1,2, ^1,1, ^1,0; ^2,3; ^2,2, ^2,l) 

One of the main obstacles is that the original equations of neural dynamics are not written in 
any of the canonical forms for which the reconstruction algorithms are available. The question, 
therefore, is if there exists an invertible coordinate transformation such that the model equations 
can be rendered canonic. Below we demonstrate that this is generally not the case if the transfor- 
mation is parameter-independent. This is formally stated in Section 3.1. However, if we allow our 
transformation to be both parameter and time-dependent, a relevant class of models with polyno- 
mial right-hand sides can be transformed into one of the canonic forms. This is demonstrated in 
Section 3.2. 



3.1. Parameter-independent time-invariant transformations 

Let us consider a class of systems that can be described by (13.41) . Clearly this system is not in a 
canonical adaptive observer form because A{\) depends on the unknown parameter A explicitly. 
The question, however, is if there exists a differentiable coordinate transformation 

z = Zi = Xi 

such that in the new coordinates the equations of system ( 13.41) satisfy one of the canonic descrip- 
tions. We show that the answer to this question is negative, and it follows from the following 
slightly more general statement 

Theorem 6. The system 

p 

X = f{x) + ^eiqi{x) 



i=l 



X E W\y e R,qi 



y = Kx) 



with 



fix) 



/ 1 1 ... 1 \ 
... 



X, 



h{x) = ( 1 ... ) 



X 



(3.6) 



(3.7) 



\ ... / 

cannot be transformed by diffeomorphic change of coordinates, z 



[xj, into 



Aiz + My) + Yl ^i^i(y)^ y = '^^^ 



■1=1 



2 e M", y G M, : M ^ W (3.8) 

with Ai^Ci in canonical observer form ( 12.751) , if either (i) n > 2 or (ii) there exists i G 
{1, ...,p}, j G {2, ...,n} such that dq-Jdxj ^ 0. 

The proof of Theorem [6] and other results are provided in the Appendix. 
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3.2. Parameter-dependent and time-varying transformations 

Let us now consider the case in which the transformation z = A, 9, t) is allowed to depend 
on unknown parameters and time. As we show below, this class of transformations is much more 
flexible. In principle it allows us to solve the problem of partial state and parameter reconstruc- 
tion for an important class of oscillators with polynomial right-hand side and time-invariant time 
constants. 

We start by searching for a transformation $: 

$: g = r(A)x, |r(A)|7^o 

such that 

T{\)A{X)T~\X) ={ll) (3-9) 

where the matrix A{\) is defined as in (13.51) . It is easy to see that the transformation satisfying this 
constraint exists, and it is determined by 

T(A)=(^ J)- (3.10) 

According to (13.91) . (13.101) and (13.51) equations of (13.41) in the coordinates q can be written as 

g = Aig + V^(gi)r7(e,A), (3.11) 



where 



' '^^'^^^ ^ gi 1 ; ' (3.12) 



rj{6, A) = (6'i_3, 6'i_2, ^1,1 — A, ^i^o, ^^1,3 + ^2,37 -^^1,2 + ^'2,27 -^^1,1 + ^2,17 -^^1, 



0) 



Remark 7. Notice that 

• availability of the parameter vector t] in l\3.11\) , A3.12\) , expressed as a function of 9, A, 
implies the availability of 6, A z/6'1 7^ 0. Indeed, in this case the value of X = rjs/rj^ and 
the values of all 9ij are uniquely defined by r]i; 

• condition ^10 7^ is sufficient for reconstructing the values ofx2 provided that q and t] are 
available; indeed in this case x = T~^(A)g 

As follows from Remark|7]the problem of state and parameter reconstruction of (13.41) from mea- 
sured data xi{t) amounts to solving the problem of state and parameter reconstruction of (13.1 II) . 
In order to solve this problem we shall employ yet another coordinate transform: 



zi = gi 

Z2 = q2 + C{t)r] 



T,.^ (3.13) 
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in which the functions ({t) are some differentiable functions of time. Coordinate transformation 
(13.131 ) is clearly time-dependent. The role of this additional transformation is to transform the 
equations of system (13.1 II) into the form for which a solution already exists. 

Definitions of these functions, specific estimation algorithms and their convergence properties 
are discussed in detail in the next section. 

3.3. Observers for transformed equations 
3.3.1. Bastin-Gevers Adaptive Observer 

Proceeding from (13.1 II) . (13.121) and applying a second change of coordinates given by 



where / G M<o and A; G M are some design parameters, we obtain the canonical form (12.71) 
presented in [|3l 



System (13.141 ) now is in the Bastin-Gevers adaptive observer canonical form. Notice that the 
parameter vector ?7(A, 9) remains unchanged and recall Remark|71 Let us proceed to the observer 
construction following the steps described in (12.81) - (12.131) . 

We start by introducing an auxiliary filter of which the general form is given by (I2.10|) . Ac- 
cording to (13.141) the auxiliary filter is defined as follows: 




T]{9, A) = {013, Oi2, Oil — A, 6^10, 6*23 + A6'i3, 6^22 + A^i2, 021 + A6'ii, A6'io)^ 




(3.14) 



Vl = fvi + 

V2 = fv2 + iy^ 

V3 = fv3 + iy 

V4 = fvi + i 

V5 = fv5 + 

v& = fve + 

V7 = fvj + \y 



(3.15) 
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Hence in accordance with (12.111) the regressor vector (p{t) is written as 

(fi = kvi + 

ip2 = kv2 + 

(P3 = kvs + y 

ip4^ = kV4^ + 1 

ip5 = kvi 

ipa = kv2 

V97 = kvi, 

V98 = kVi 

and the observer equations are as follows: 



(3.16) 



17 = Tif{xi - xi), 
V" = (wi, W2, t's, ^4, ^'s, ^'e, vj, vs) 



(3.17) 



Taking (13.151 ) - (13.171 ). and (12.131 ) into account we obtain the following equations governing 
the dynamics of the estimation error, (x*, f/)^ 

""^^ ( 0' / )'^*+ ( '^0'' ) (3.18) 
7] = —Tipxl 

The auxiliary filter (13.151) acts here as an inherent component of a time-varying coordinate trans- 
formation rendering the error dynamics into (12.131) . This coordinate transformation is similar to 
that defined by (13.131) . provided that z, i] in (13.131) are replaced by estimation errors x*, f]. 

Let us now explore asymptotic properties of the observer. First we notice that f4(t), v^{t) both 
converge to constant values exponentially fast as t —t- 00. In fact, 

lim Vi{t) = - J, lim vs^t) = 

t^oo K t^oo JK 

Thus accordingly ^Pi{t), (psit) both tend to constant values as t — )■ 00: 

lim ip4{t) = 0, lim (psit) = 

t— >oo t—>oo f 

The latter fact implies that the persistency of excitation requirement is necessarily violated for 
regressor (13.161 ). Indeed, condition (12.141) does not hold if one of the components of (p{t) is ex- 
ponentially converging to zero. The question therefore, is if this approach can be used at all to 
construct asymptotically converging estimators of state and parameters of (|3.14l) . The answer to 
this question is provided in the corollary below 
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Corollary 8. Consider the function 

If it is globally bounded and persistently exciting, (\2.14\l , then the following holds along the solu- 
tions of Bl5\i - Bl7\i : 

lim x{t) — x{t) = 0, lim r)j(t) — T]i, i ^ 4, 

t—^oo t—^oo 

and the convergence is exponential. 

Remark 9. Corollary |5] demonstrates that despite the original result of /. e. Theorem |?1 does 
not apply to system ([3.14\l directly one can still construct a reduced order observer for this system. 
This reduced observer guarantees partial reconstruction of unmeasured parameters, and this re- 
construction is exponentially fast. To recover the true values of unknown parameters one needs to 
solve the following system 

Vl = ^13 
V2 = dl2 
Vi = 0u- \ 

V5 = O23 + A6'i3 
Ve = 022 + A6'i2 
Vr = 021 + A6'ii 
r]s = \9io 

for 6i, A taking the values offji as the estimates ofrji. Solution to this system may not be unique, 
hence the reconstruction is generally possible only up to a certain scaling factor 

Simulation results for this observer are presented in Section 5. 
3.3.2. Marino-Tomei Observer 

Let us define the vector-function C(t) in (13.131) as follows: 

Ci = -Ki + ktpiAii) - ^2,i(gi), k e ]R>o 

In this case we have 

8 / ^ \ * 

i2 = ^ ^2,i(gi)?7i + k l-^Qrji + i)iA{qi)r]i | - ^ ^/'2,i(gi)?7i 
i=i \ i=i J i=i 

8 

= k^{-Ci + iJiMi))vi 

i=l 
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Hence, taking equality (13.131) into account and expressing q2 as q2 = Z2 — C'^{t)7]{6, A) we obtain 

8 

Zi = Z2 + ^{-Ci + AAll))Vi 
i=l 

8 

Z2 = k ^(-C + iJi,i{Qi))Vi 



(3.19) 



i=l 



Notice that iJi^iiqi) = i^ifiiQi) = 1. hence — C4 + ipi,A{qi) and — + i^ifiiQi) converge to some 
constants in M exponentially fast as t — )• 00. Moreover, the sum —^4 + ^^1.4(^1) is converging to 
zero, and the sum —C^ + ^1.8(91) is converging to — 1/A; as t — 00. Taking these facts into account 
we can conclude that system (13.191) can be rewritten in the following (reduced) form 

z = Aiz + b(f)'^{zi,t)v{6, A) + be{t), 



A, 



(3.20) 



1 


/ -Ci + V^i,i(^i) \ 
-C2 + ^^1,2(^1) 
-C3 + ^^1,3(^1) 
-C8 + ^lA^i) 

-C5+^l,5{zi) 

\ -Ct + ^1,7(^1) / 
v{6, A) = (6'i^3, 9i^2, 01,1, ^1,05 -^^1,3 

where e{t) is an exponentially decaying term. 

System (13.201 ) is clearly in the adaptive canonic observer form. Hence it admits the following 
adaptive observer 

Aiz + L{z - z) +b(l)^{zi,t)v 
-h ~ 



^2,35 -^^1,2 + ^2,2, ^01,1 + 02,1)^ 



Z 

L 



li = k + l,l2 = k 



(3.21) 



-h 

Vi = -7(^1 - zi)(f)i{zi,t), 7 G ]R>o 
of which the asymptotic properties are specified in the following Theorem 

Theorem 10. Let us suppose that system ( I3.20P be given and its solutions are defined for all t. 
Then, for all initial conditions, solutions of the combined system ^3.201) . ^3.2iP exist for all t and 

lim z{t) - z{t) = 

Furthermore, if the function (f){zi, t) is persistently exciting and z{t) is bounded then 

lim v{t) - v{e,X) = 0, 

and the dynamics of z — z,v — v are exponentially stable in the sense of Lyapunov. 
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The proof of Theorem [TO] is provided in the Appendix. 

Remark 11. Similar to Corollary \8\ for Bastin-Gevers observer, Theorem U0\ provides us with a 
computational scheme that, subject to that (j){zi,t) is persistently exciting, can be used to estimate 
the values of the modified vector of uncertain parameters v{6, A). The question, however, is that if 
the values of 9, A can always be restored from v{9, A). In general, the answer to this question is 
negative. Indeed, according to di.iil) we have 



As follows from {[3.22\) one can easily recover the values o/6'i 3, 6^1 2, and 9i^i. However, recovering 
the values of remaining parameters explicitly from the estimates ofv{9, A) is possible only up to 
a certain scaling parameter. Indeed, if the number of unknowns in (13.22\) exceeds the number of 
equations by one. 

Remark 12. Notice that in the relevant special cases, when the value of either 6^2,3, 6^2, 2, or 6*2, 1 is 
zero, such reconstruction is obviously possible. Let us suppose that 6*2,3 = 0. Hence the value ofX 
can be expressed from A3.22\i as 



and thus the rest of parameters can be reconstructed as well. Due to the presence of division in 
(13.23\) , this scheme may be sensitive to persistent perturbations when = A6'i 3 is small. 

So far we considered special cases of (11.31) in which the time constants of unmeasured vari- 
ables were unknown yet constant and parametrization of the right-hand side was linear. As we 
mentioned in Remark [T2l even for this simpler class of systems solving the problem of parameter 
reconstruction may not be a straightforward operation. For example, if there are cubic, quadratic 
and linear terms in the second equation of (13.41) then recovering all parameters of (13.41) by observer 
(13.211) may not be possible. Nonlinear parametrization, time-varying time constants and nonlin- 
ear coupling between equations in the right-hand side of (11.31 ) make the reconstruction problem 
even more complicated. Even though there are results that partially address the issue of nonlinear 
parametrization, see e.g. |[34ll , [[331 , |[6l , |[32l , [[T8l , the estimation problem for systems with general 
nonlinear parametrization is still an open issue. 

In the next section we show that for a large subclass of (11.31) there always exists an observer 
that solves the problem of state and parameter reconstruction from the measurements of V. More- 
over the structure of this observer does not depend significantly on specific equations describing 




(3.22) 



A 



(3.23) 
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dynamics of the observed system. For this reason, and similarly to [[TTll . we refer to this class of 
observers as universal adaptive observers. 



4. Universal adaptive observers for conductance-based models 

The ideas of universal adaptive observers for systems with nonlinearly parameterized uncertainty 
was introduced in a series of works Il35]| . [|36]| devoted to the study of convergence to unstable 
invariant sets. Here we provide a review of these results and discuss how they can be applied to 
the problem of state and parameter reconstruction of (11.31 ). 
The following class of models is considered in [|36l : 

xo =9^(f)oixo,Po,t) + Y,'^=iCi{xo,qi,t)xi + Co{xo,qo,t)+^o{t)+u{t) 
ii = -(3i{xo,Ti,t) xi + 9f(pi{xo,pi,t) 

Xi = -l3i{xo,Ti,t) Xi + 6f(j)i{xo,pi,t) + ^i{t), 

^ Xn = -(3niXo,rn,t) Xn + 0[l^(l)i{Xo,Pn,t) + ^n{t), 



are continuous and known functions, u : R>o — M, m G C° is a known function of time modelling 
the control input, and : R>o — M, G C° are functions that are unknown, yet bounded. The 
functions ^j(t) represent unmodeled dynamics, external perturbations, residuals due to the coarse- 
graining procedures at the stage of reduction [8], etc. 

Variable y in system (14.11) is the output, and the variables Xi,i>l are the components of state 
X, that are not available for direct observation. Vectors 6i G consist of linear parameters of 
uncertainties in the right-hand side of the i-th equation in (|4.1I) . Parameters r^GlR, « = {!,. ..,«} 
are the unknown parameters of time-varying relaxation rates, /3i(xo, Tj, t), of the state variables Xj, 
and vectors pi G R™% G W\ consist of the nonlinear parameters of the uncertainties. The 
functions q(xo, qi, t) are supposed to be bounded. 

Notice that system (14.11) is almost as general as (11.31) . The only difference is that variables Xi, 
i>l enter the first equation of (14.11) as 



y = (1, 0, . . . , 0)x = Xo, Xiito) = Xi^o G M, 

X = Col(Xo, Xi,...,Xn), 9i = C0l(6'i_i, . . . , Oi^di), 



where 



0i :R X R™' X R>o M*, 0i G C°, di, G N, i = {0, . . . , n} 

A :R X R X M>o M>o, A G C°, z = {1, . . . , n} 

d :R X R'' X R>o ^ R, q G (3°, G N, i = {0,...,n} 



n 




i=l 
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whereas the corresponding variables in system (11.31 ) enter the first equation in a slightly more 
general way 

This difference, however, is not critical for the observers presented in ||36l can be adjusted to deal 
with this more general case as well. 

For notational convenience we denote: 

X = Col(po,gO,n,Pl,gi • • ■^Tn.PnAn), 

n 

s = dim (A) = n + ^^(mj + r^). 

Symbols and Qx, respectively, denote domains of admissible values for 6 and A. 

The system state x = col(xo, xi, - ■ ■ , x„) is not measured; only the values of the input u(t) 
and the output y{t) = Xo{t),t > to in (14.11 ) are accessible over any time interval [to, t] that belongs 
to the history of the system. The actual values of parameters 6, A are assumed to be unknown 
a-priori. We assume however, that they belong to a set, e.g. a hypercube, with known bounds: 

Instead of imposing the traditional requirement of asymptotic estimation of the unknown pa- 
rameters with arbitrarily small error we relax our demands to estimating the values of state and 
parameters of (14.11 ) up to a certain tolerance. This is because we allow unmodeled dynamics, ^i{t), 
in the right-hand side of (14.11 ). As a result of such a practically important addition there may exist a 
set of systems of which the solutions are relatively close to the measured data yet their parameters 
could be different. Instead of just one value of unknown parameter vectors 0, A we therefore have 
to deal with a set of 0, A corresponding to the solutions of (14.11 ) that over time are sufficiently 
close. This set of model parameters is referred to as an equivalence class of (14.11 ). 

Similarly to canonical observer schemes [|23l , O, [|25l the method presented in [|36l relies on 
the ability to evaluate the integrals 

/^.(t,r,,Pi) ^ f e-~ ^rP^iMx),n,x)dx^^^^^^^)^p^^^)dT (4.2) 

J to 

at a given time t and for the given values of Tj, pi within a given accuracy. In classical adaptive 
observer schemes, the values of /3j(xo, Tj, t) are constant. This allows us to transform the original 
equations by a (possibly parameter-dependent) non-singular linear coordinate transformation, $ : 
X (-7- z, xi = zi, into an equivalent form in which the values of all time constants are known. In 
the new coordinates the variables Z2, ■ ■ ■ ,Zn can be estimated by integrals (14.21) in which the values 
of /3i(xo, Tj, t) are constant and known. This is usually done by using auxiliary linear filters. In our 
case, the values of A(xo, Tj, t) are not constant and are unknown due to the presence of r^. Yet if 



23 



D. Fairhurst et al. 



Observers for Canonic Models 



the values of rj would be known we could still estimate the values of integrals (14.21 ) as follows 

/ e-^-'^'(^«(^)'^-^)'^%(xo(r),p„r)c/r ~ (4.3) 

J to 
Jt~T 

where T E M>o is sufficiently large and t >T + tQ. 

Alternatively, if 0j(xo(t),Pi, t), /3j(xo(t), Tj, t) are periodic with rationally - dependent periods 
and satisfy the Dini condition in t, integrals (14.21 ) can be estimated invoking a Fourier expansion. 
Notice that for continuous and Lipschitz in pi functions Tj, pj) the coefficients of their Fourier 
expansion remain continuous and Lipschitz with respect to pi. 

In the next sections we present the general structure of the observer for (14.11) and provide a list 
of its asymptotic properties. 

4.1. Observer definition and assumptions 

Consider the following function (^(xq, A, t) : R x x R>o -^R'^,d = X;r=o ^i' 

<^(xo, A,t) = (0o(xo,Po,^),ci(xo,gi,t)/ii(t,ri,pi), . . . 

T (4.4) 

• • • ,Cn{Xo,qn,t)fin{t,Tn,Pn)) 

The function (p{xo, A, t) is a concatenation of 0o(') and integrals (14.21) . We assume that the values 
of (p{xo, A, t) can be efficiently estimated for all Xq, A, t > up to a small mismatch. In other 
words, we suppose that there exists a function ^(xq, A, t) such that the following property holds: 

||c^(xo,A,t) -cp(xo,A,t)|| < A^, A^eM>o, (4.5) 

where values of ^(xq, A, t) are efficiently computable for all xq. A, t (see e.g. (14.31) for an example 
of such approximations), and is sufficiently small. 

If parameters Ti,pi, and qt in the right-hand side of (14.11) would be known and Cj(xo, qi,t) = 1, 
(3i{xo, Ti, t) = Ti, then the function ip{xo, A, t) could be estimated by {(j)o{xo, t),r]i, . . . , rfn) where 
rji are the solutions of the following auxiliary system (filter) 

m = -^iVi + Mxo,Ph t) (4.6) 

with zero initial conditions. Systems like (14.61 ) are inherent components of standard adaptive ob- 
servers [[T71 , ||3l , [|23l . In our case we suppose that the values of t^, qi, pi are not know a-priori and 
that Cj(xo, qi,t), l3i{xo, Tj, t) are not constant. Therefore, we replace rji with their approximations, 
e.g. as in (14.31) : 

(fi{xo,X,t) = (0o(xo,po,^),ci(xo,gi,t)/ii(t,ri,pi), . . . 

■ ■ ■ : C-i^i^XQ, qn,t^ fin(ty Tfi^ Pn)) 
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For periodic 0j(xo(t),Pi, t), /3j(xo(t), r,, t) a Fourier expansion can be employed to define (p{xo, A, t). 
The value of in (14.51 ) stands for the accuracy of approximation, and as a rule of thumb the more 
computational resources are devoted to approximate ifi{xo, A, t) the smaller is the value of A^. 

With regard to the functions ^i{t) in (14.11) we suppose that an upper bound, A^, of the following 
sum is available: 

n ^ 

J2 -\Mr)\\oo,ito,o.] + IKo(r)||oo,[to,oo] < A^, A^ G M>o. (4.7) 

• 1 i 

Denoting co{xo,qo,t) = Co{xo,X,t), for notational convenience, we can now define the ob- 
server as 

xo = -a(xo -xq) + 6 (p{xo, A, t) + co(xo, A, t) + u{t) 
= -leixo - Xo)ip{xo, A, t), 70, a G M>o 
Xi = -(3i{xo,fi,t)xi + 9j(l)i{xo,Pi,t), i = {l,...,n}, (4.9) 

where 

e = co\ieo,ei,--- ,0n) 

is the vector of estimates of 6. The components of vector A = col(po, ^0, ^-1,^1, gi, . . . , fn,Pn, Qn) = 
col(Ai, . . . , As), with s = dim (A), evolve according to the following equations 




7 ■ ■ e ■ {xij - X2J - xij [xlj + xlj)) 

7 ■ CUj ■ e ■ [xij + X2,j - X2,j {xj j + xlj)) 

\ _|_ -^j.max^-^j.min / _|_ 1 

a{\\xo - Xolle), 



(4.10) 



j = {l,...,s}, xl^ito) + xl^ito) = 1, (4.11) 

where a{-) : M — t- ]R>o is abounded continuous function, i.e. a(v) < S E M>o, and \o'{v)\ < \v\ 
for all t; G M. We set Uj G M>o and let Uj be rationally-independent: 

'^tOjkj ^0, \/ kj eZ. (4.12) 



In order to proceed further we will need the notions of X-uniform persistency of excitation [|2T 
and nonlinear persistency of excitation [[6l: 



Definition 13 (A-uniform persistency of excitation). Let (p : R>o x V ^ M"^"", V C W be a 
continuous function. We say that ifi{t, A) is X-uniformly persistently exciting (X-uPE) if there exist 
fi G M>o. L G M>o such that for each \ eT) 

rt+L 

cp{t, X)(p{t, xfdr >fiim> to. (4.13) 
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In contrast to conventional definitions, the present notion requires that the lower bound for the 
integral J^^^ cp{t, X)^{t, X)'^dT in (14.131) does not vanish for all X G V, and is separated away 
from zero. We need this property in order to determine the linear parts, 9i, of the parametric 
uncertainties in model (14.11) . 

To reconstruct the nonlinear part of the uncertainties. A, we will require that (fi{xo, X, t) is non- 
linearly persistently exciting in A. Here we adopt the definition of nonlinear persistent excitation 
from [|6l with a minor modification. The modification is needed to account for a possibility that 

ifixo, A, t) = ifixo, A', t), A ^ A', t G M, 

which is the case, for example if ifi{xo, A, t) is periodic in A. The modified notion is presented in 
Definition [141 below. 



Definition 14 (Nonlinear persistency of excitation). The function 

(p{xo, A, t) is nonlinearly persistently exciting if there exist L, 13 E ]R>o such that for all A, A' G Qx 
and t G M there exists t' G \t — L,t] ensuring that the following inequality holds 

||(^(a;o,A,t)-(^(xo,A',t')ll dist(£'(A),A'), (4.14) 

^(A) = {A' G ^^aI ^{xq, a', t) = cp{xo, A, t) V t G M} (4.15) 

The symbol S{X) denotes the equivalence class for A, and dist(£'(A), A') in (14.141 ) substitutes 
the Euclidian norm in the [61 original definition. The nonlinear persistency of excitation condition 
( 14.141) is very similar to its linear counterpart (14.131 ). In fact (14.131 ) can be written in the form of 
inequality (14.141) . cf. fT7\. For further discussion of these notions, see ETIl . 



4.2. Asymptotic properties of the observer 

The main results of this section are provided in Theorems \T5\ and [171 Theorem \T5\ establishes 
conditions for state boundedness of the observer, and states its general asymptotic properties. The- 
orem [l71 specifies a set of conditions for the possibility of asymptotic reconstruction of 9i, Ti, and 
Pi, up to their equivalence classes and small mismatch due to errors. 

Proofs of Theorems [TSl [TTl and other auxiliary results can be found in [|36ll . 



Theorem 15 (Boundedness). Let system M.Ik ( 14.^1) - ([4.10^ be given. Assume that function 
(p{xo{t), A, t) is \-uniformly persistently exciting, and the functions ip{xo(t), A, t), co(xo(t), A, t) 
are Lipschitz in X: 

mxo{t),\t) - (^(xo(t), A',t)|| < D\\X->^1 
||co(xo(t),A,t) -co(xo(t),A',t)|| < DJA- A'||. 

Then there exist numbers £ > 0, 7* > such that for all 7 G (0, 7*].' 
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1 ) trajectories of the closed loop system ([4.8\l - A4.10\) are bounded and 

lim ||xo(t)-Xo(t)||, = 0; (4.17) 



2) there exists X* G Q\, k G ]R>o such that 

lim \{t) = X* 

t—>oo 

limsup ll^(t) - e\\ < K{{D\\e\\ + D,)\\X* - A|| + 2A). 

t—^oo 

A = ||0||A^ + A^ (4.19) 



(4.18) 



Remark 16. TheoremfT?] assures that the estimates 0{t), X{t) asymptotically converge to a neigh- 
borhood of the actual values 0, X. It does not specify, however, how close these estimates are to 
the true values of 6, X. 

The next result states that if the values of A,„ and At: 



\(p{xo,X.t) - <pixo,X,t)\\ < A 



1 

^-\\U^)\\oo,[to,oo] + ||^o(r)||oo,[to,oo] < 

in (14.51) . (14.71) are small, e.g. ^(xo(t), A, t) approximates ip{xo{t), X,t) with sufficiently high 
accuracy and the unmodeled dynamic is negligible, the estimates 9(t), X{t) will converge to small 
neighborhoods of the equivalence classes of 0, X. The sizes of these neighborhoods are shown to 
be bounded from above by monotone functions A: 

A = ||0||A^ + A^ 

vanishing at zero. Formally this result is stated in Theorem [17] below 

Theorem 17 (Convergence). Let the assumptions of Theorem\T5\hold, assume that <^o(2^0) \ t) G 
C^, the derivative d(^Q{xo{t), X, t)/dt is globally bounded, and A = ||^|| A^^ + A^ is small. Then 
there exist numbers e > 0, 7* > such that for all 7j G (0, 7*) 

1) 

limsup \\0{t) - 0\\ = 0{VA) + C(A) 

2) in case 0^ (p{xo{t), A, t) + co(xo(t), A, t) is nonlinearly persistently exciting with respect to 
X, then the estimates A(t) converge into a small vicinity of £{X): 

limsup dist(A(t),^(A)) = O(v^) + C(A) (4.20) 
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5. Examples 

5.1. Parameter estimation of the 2D Hindmarsh-Rose model with Bastin- 
Gevers observer 

The canonical form (|3.14l) and the observer presented in Section 3.3 were built in MATLAB and 
using the differential equation solver ode45, numerical results were obtained. Figure [3] shows the 
parameter convergence of each r]i{t), i ^ A. The various parameter values were set as follows: for 
the neuron 

A = 2.027,^13 = -10.4,^12 = -4.35,^11 = 6.65, ^lo = 0.9125,^22 = -32.45, = -32.15 
and for the observer dynamics 

A; = 1, ci = 1, r = diag(l, 1, 1, 1, 1, 1, 1, 1), F = -1 



5.2. Parameter estimation of the 2D Hindmarsh-Rose model with Mario- 
Tomei observer 

In addition to simulating the Bastin-Gever observer for the Hindmarsh-Rose model we also simu- 
lated observer (I3.21I) derived within the framework of the approach presented in ||23l . In this case 
we set 

^1,3 = -1, ^1,2 = 3, ^1,3 = 0, ^1,0 = 1.5, 02,3 = 0, ^2,2 = "5, ^2,1 = 0, A = -1 

There is no particular reasoning behind our choice of parameters in both this and previous case 
apart from that these parameters must induce persistent oscillatory dynamics of the solutions of 
(I3.3|) . Parameters of the observer were chosen as follows: 

h = k = I, 7 = 1, k = 1. 
Simulation results for this system are shown in Figure SI 

5.3. Parameter estimation of the Morris-Lecar model 

Let us now turn to a more realistic class of equations, i.e. conductance-based models. In particular, 
we consider the Morris-Lecar model, ||28l : 

V = ^ {-9c.m^{V){V - Eca) - gKw{V - Ek) - - ^o)) + / 

1 , ^oo(V^) ^^-^^ 
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2r 8r 10 




(g) fis{t)w. t 

Figure 3: Simulation results for Bastin-Gevers observer. Each T)i{t),i ^ A are shown with the true 
values indicated with a broken line. The periodic time of the spikes is circa 10 seconds while the 
total time simulated is 5 x 10^ seconds, thus the figure spans thousands of spike cycles. The visual 
effect of these extreme time scales is seen in the figure; hunting oscillations within the observer 
are seen as thickening of the graphs of certain estimated parameters. 
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(b) V2it) V. t 



2 4 6 
(d) Vi{t) V. t 



8 x10' 



2 4 6 
(e) v^it) V. t 



(C) V3{t) V. t 



1 -4 

x10^ 



2 4 6 
if) ve{t) V. t 



x10 



2 4 6 8 x 10' 
(g) Vrit) V. t 

Figure 4: Simulation results for Marino-Tomei observer. Each Vi{t) are shown with the true values 
indicated with a broken line. The periodic time of the spikes is circa 10 seconds while the total 
time simulated is 10 x 10'^ seconds, thus the figure spans thousands of spike cycles. The visual 
effect of these extreme time scales is seen in the figure; hunting oscillations within the observer 
are seen as thickening of the graphs of certain estimated parameters. 
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where 

moo{V) 
t{V) 

System (15.11) is a reduction of the standard 4-dimensional Hodgkin-Huxley equations, and is one 
of the simplest models describing the dynamics of evoked membrane potential and, at the same 
time, claiming biological plausibility. 

Parameters gcn, gx, and gi stand for the maximal conductances of the calcium, potassium and 
leakage currents respectively; C is the membrane capacitance; Vi, V2, V3, V4 are the parameters of 
the gating variables; Tq is the parameter regulating the time scale of ionic currents; Ec^, and are 
the Nernst potentials of the calcium and potassium currents, and E^ is the rest potential. Variable 
/ models an external stimulation current. In this example the value of / was set to / = 10. 

The total number of parameters in system (15.11) is 12, excluding the stimulation current /. Some 
of these parameters, however, are already available or can be considered typical. For example the 
values of the Nernst potentials for calcium and potassium channels, Eca, Ek, are known and 
usually are set as follows .Eca = 100, E^ = —70. The value of the rest potential, Vq, can be 
estimated from the cell explicitly. Here we set Vq = —50. Parameters Vi, V2 characterize the 
steady-state response curve of the activation gates corresponding to the calcium channels, and V3, 
V4 are the parameters of the potassium channels. In the simulations we set these parameters to 
standard values as e.g. in [[Ml: Vi = -1, V2 = 15, ^3 = 10, and V4 = 29. 

The values of parameters, gca., giK, dL, and Tq, however, may vary substantially from one cell 
to another. For example, the values of gcs,, gx, Ql depend on the density of ion channels in a 
patch of the membrane; the value of Tq is dependent on temperature. Hence, in order to model the 
dynamics of individual cells, we need to be able to recover these values from data. 

As before, we suppose that the values of V over time are available for direct observation, and 
the values of w are not measured. System (15.11) has no linear time-invariant part, and the dynamics 
of w are governed by a nonlinear differential equation with the time-varying relaxation factor, 
l/riy). Therefore, observers presented in Section 3 may not be applied explicitly to this system. 
This does not imply, however, that parameters gcg,, cjk, gi, and Tq cannot be recovered from the 
measurements of V . In fact, as we show below, one can successfully reconstruct these parameters 
by using observers defined in Section 4. 

For the sake of notational consistency we denote Xq = V, Xi = w, and without loss of gener- 
ality suppose that (7 = 1. Hence system (15.11 ) can now be rewritten as follows: 

io = do,m^{xo){xo - ^ca) + 6'o,ia;i(a;o - Ek) + 9o,3{xo - Vq) + / 

1 (5.2) 

xi = -/3i(a;o, A)xi + -r0i(xo) 



0.5 1 + tanh 



0.5 I 1 + tanh 
1 



V -Vi 
V-V3 



Tr 



0" 



cosh 



V-V3 
2V4, 
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where 



/3i{xo, A) = Y cosh ( ^° ) , A = To 



A \ 2V4 



'A 



Noticing that /3i(xo, A) is separated away from zero for all bounded xq and positive A we substitute 
variable Xx in (15.21 ) with its estimation 

The larger the value of T the higher the accuracy of estimation for large t. After this substitution 
system (15.21 ) reduces to just only one equation 

io = 6'o,i0o,i(a;o) + 6^0,200,2 (a^o, A, t) + 6^0,300,3 (a^o) + I + ^o{t) (5.3) 

where 

00,1 (a;o) = "^oo(a;o)(a;o - -^ca) 
0o,2(a;o, A, t) = (xq - ^x)x(A, t) 
00,3 (a^o) = {xo - Vo) 

and term ,^o(^) is bounded. 

Equation (|5.3I) is a special case of (14.11) . and hence we can apply the results of Section 4 to 
construct an observer for asymptotic estimation of the values of 6*0,1, 6*0,2, 6*0,3, and A. In accordance 
with (14. 8|) - (|4.10l) we obtain the following observer equations 



= -a{xo - xo) + 6'i0o,i(a;o) + 6'20o,2(a;o, A, t) + 6'30o,3(a;o) + / 

01 = -le{xo - a;o)0o,i(xo) 

02 = -le{xo - a;o)0o,2(a;o, A, t) 
^3 = -leixo - a;o)0o,3(xo) 

A = 3 + 

£1,1 = le - £2,1 - 1 + 

£2,1 = le (xi,i + X2,i - X2,i [xl i + X2,i)) 
e = (t(||xo - Xolle) 



(5.4) 



(5.5) 



Parameters of the observer were set as follows: a = l,e = 0.001, 7 = 0.01, 70 = 0.05. 

According to Theorems [121 [17] observer (15.41) . (15.51) should ensure successful reconstruction 
of the model parameters provided that the regressor is persistently exciting. This requirement is 
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2 4 6 x104,? 2 4 6 xlO^,? 



Figure 5: Trajectories of the estimates of parameters 6i, 62, 6^ and A as functions of time. The 
periodic time of the spikes is circa 10 seconds while the total time simulated is 8 x 10"^ seconds, 
thus the figure spans thousands of spike cycles. The visual effect of these extreme time scales is 
seen in the figure; hunting oscillations within the observer are seen as thickening of the graphs of 
certain estimated parameters. 

satisfied for model (15.21) generating periodic solutions. We simulated system (15.21) . (15.41) . (15.51) over 
a wide range of initial conditions. Figure [5] shows an example of typical behavior of the observer 
over time. As we can see from this figure all estimates converge to small neighborhoods of true 
values of the parameters. 



6. Conclusion 

In this article we have reviewed and explored observer-based approaches to the problem of state 
and parameter reconstruction for classes of typical models of neural oscillators. The estimation 
procedure in this approach is defined as a system of ordinary differential equations of which the 
right-hand side does not depend explicitly on the unmeasured variables. The solution of this system 
(or functions of the solutions) should asymptotically converge to small neighbourhoods of the 
actual values of the variables to be estimated. Until recently, due to nonlinear dependence of 
the vector-fields of the models on unknown parameters and also due to uncertainties in the time 
scales of hidden variables, observer-based approach to solving the problem of state and parameter 
estimation of neural oscillators was a relatively unexplored territory. Here we demonstrate that 
despite these obvious difficulties the approach can be successfully applied to a wide range of 
models. 

Two different strategies to observer design have been studied in the paper. The first strategy is 
based on the availability of canonical representations of the original system. Success of this strat- 
egy is obviously determined by wether one can find a suitable coordinate transformation such that 
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the equations of the original model can be transformed into the canonical adaptive observer form. 
Because a coordinate transformation is required, different classes of models are likely to lead to 
different observers. The second strategy is based on the ideas and approaches of universal adaptive 
regulation ifTTI . non-uniform convergence and non-uniform small-gain theorems |[35l . Il36ll . The 
structure of observers obtained as a result of this design strategy does not change much from one 
model to another. The main difference between these design strategies is in the convergence rates: 
exponential for the first and asymptotical for the second. As long as mere overall convergence time 
is accounted for there is no big difference whether the convergence itself is exponential or not. Yet, 
the fact that it can be made exponential with known rates of convergence allows us to derive the 
a-priori estimates of the amount of time needed to achieve a certain given accuracy of estimation. 

We have shown that for linearly parameterized models such as the FitzHugh-Nagumo and 
Hindmarsh-Rose oscillators one can develop an observer for state and parameter estimation of 
which the convergence rate is exponential. For the nonlinearly parameterized and more realistic 
models such as the Morris-Lecar and Hodgkin-Huxley equations we presented an observer of 
which the convergence is asymptotic. In both cases the rate of convergence depends on the degree 
of excitation in the measured data. In the case of linearly parameterized systems this excitation 
can be measured by the minimal eigenvalue of a certain matrix constructed explicitly from the data 
and the model. For the nonlinearly parameterized systems the degree of excitation is defined by a 
more complex expression, (14.141 ). In principle, one can ensure arbitrarily fast convergence of the 
estimator provided that the excitation is sufficiently high. This property motivates the development 
of measurements protocols that are most consistent with a range of models that will be fitted to the 
collected data. In fact, in order to achieve higher computational effectiveness, one shall aim to 
produce data of which the excitation is higher for the given range of models. 

One question remains unexplored though - the actual amount of elementary computational 
operations required to realize these two observer schemes. This number depends substantially 
on the required accuracy of estimation. We aim to answer this important question in future case 
studies. 
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7. Appendix 

Proof of Theorem |6| 

Proof. The proof is straightforward. Indeed, for the observability test we have 
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(7.2) 



V ^-^h{x) j 

It follows that observability is lost for n > 2. This proves condition (i). In order to demonstrate 
condition (ii) we use theorem [5] as follows. From equation (12.211 ) we have 



1 
1 



9\ 

92 



(7.3) 



giving 



(7.4) 
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Theorem [5] condition (i) will be satisfied since the system is linear and from theorem [5] part (ii) we 
have 

[qi,a(fjg] = [qi, g] (7.5) 

(7.6) 

(7.7) 
(7.8) 



[<li,9] 




dg 


dqi 


dx 


dx 


dqi / 


n 


dx \ 





dqi 



8X2 

we satisfy theorem [5] part (ii) if and only if dqi/ 8x2 = for all x G f/. □ 
Proof of Corollary \8\ 

Proof. The proof of the corollary is straightforward. Consider the error system given by ( 13.181) . 
Given that the function Lp(t) is bounded one can easily see that fj, x are bounded as well (consider 
e.g. the following Lyapunov candidate: V = + fiT~^f]). This implies that component 

ipi{t)fji{t) is converging to zero exponentially. 
Let us denote 

and consider the following reduced error dynamics 

r=(-;' *)^-h-(^o')h-p.W p.,) 

r] = —T(pxl 

in which e{t) stands for the term (yC4(t)f/4(t), and p = {1, 0)^. System (17.91 ) is a linear time-varying 
system of which the homogenous part is exponentially stable provided that Lp{t) is persistently 
exciting (this follows explicitly from Theorem H]). Hence, taking into account that E{t) is an expo- 
nentially converging to zero term, we can conclude that x*, fj converge to the origin too and that 
such convergence is exponential. □ 

Proof of Theorem [70] 

Proof. The proof of the theorem is standard and can be constructed from many other more general 
results (see for example ll23l . [|29ll ). Here we present just a sketch of the argument for consistency. 
According to our assumptions, matrix Ai + LCi is Hurwitz. Moreover, the transfer function 

H{s) = Ci{s-{Ai + LCi)r'h ' + ^ ' + ^ 



s"^ + {k + l)s + k (s + fc)(s + l) 

is strictly positive real. Hence, using the Kalman-Yakubovich-Popov lemma, we can conclude that 
there exists a symmetric and positive definite matrix H such that 

H{Ai + LCi) + {Al + LCifH < -Q, Hb={l, Of, (7.10) 
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where Q is a positive definite matrix. Let us now consider the following function 

1 1 f°° 

V{z,v,t) = -{z ~ zfH{z - z) + -\\v - vf-f-' + D / e\T)dT 

where the value of D is to be specified later. Clearly, the function V is well-defined for the term 
e{t) is continuous and exponentially decaying to zero as t — )■ oo. Thus the boundedness of V 
implies that p — z|| and HO — ?;|| are bounded. 
Consider the time-derivative of V: 

V = {z- zf{H{Ai + LCi) + {Al + LCifH){z - z) + {z ~ zfHbf{zi, t){v - v) 
+ {z- zfHhe -{v- {))(1, 0)(^ - 2)0(21, t) - De'^ 

Taking (17.101) into account we obtain: 

V <-{z- zfQ{z -z) + \\z - z\\\e{t)\M - Ds^ < -a\\z - z\\'^+ 
\\z- z\\\£{t)\M - De^, 



(7.11) 



where a > is the minimal eigenvalue of Q, M is fixed positive number, and D is a parameter of 
V which can be chosen arbitrarily and independently of M. Choosing D such that 



— = M, 



we ensure that 

V<--\\z~z\ 



Thus the function V is bounded from above, and given that the solution z(t) exists for all t > so 
does the solution of the combined system. The rest of the proof follows directly from Barbalatt's 
lemma and the asymptotic stability theorem for the class of skew-symmetric time- varying systems 
presented in [l27| . □ 
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